
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Standard Errors 

clear
cd ../china/data
load hs10_indlist.raw;
indN=length(hs10_indlist);

beta=zeros(indN,4);
standard_errors=zeros(indN,4);

sample_select=zeros(50,1);
for y=1:50
if (y>1 & y<5)
sample_select(y,1)=1;
else
if (y>5 & y<43)
sample_select(y,1)=1;
end
end
end
end


for w=1:50

if sample_select(w,1)==1


folder=sprintf('../china/estimation/hs10/regular/raw_%d',hs10_indlist(w,1))


cd (folder)


load imp.raw;
load data2_q_regular.mat
load beta_newq_regular.mat;
load EV_newq.mat

else

folder=sprintf('../china/estimation/hs10/biggest/raw_%d',hs10_indlist(w,1))


cd (folder)


load imp.raw;
load data2_q_biggest.mat
load beta_newq_biggest.mat;
load EV.mat

end


nExp = X;
nImp = M;
S=N*nExp*nExp;
delta=0.975;

partials=[ep switched city_switched lambda];
Fisher=zeros(4,4);

U=beta_p_sol*ep + beta_x_sol*switched + beta_c_sol*city_switched + xi_tilde_sol*lambda;

CSVF= U + delta*EV_sol;
exp_CSVF=exp(CSVF);

A=zeros(N*X*X,X);
for i=1:X
  A2=[N*X*(i-1)+1:N*X*i];
  A2=reshape(A2,N,X);
  A2=A2';
  A2=A2(:);
  A2=repmat(A2,X,1);
  A(:,i)=A2;
end
A=A(:);
A=reshape(A,X,N*X*X);

B=sum(exp(CSVF(A)),1);
B=B(:);
B_inv=1./B;

C=zeros(S,4);

for j=1:4 
  Z=exp_CSVF.*partials(:,j);
  C1=sum(Z(A),1);
  C(:,j)=C1(:);
end

M=zeros(S,16);

for x=1:4
for y=1:4

  Z=partials(:,y).*partials(:,x).*exp_CSVF;
  M1=sum(Z(A),1);
  M(:,(x-1)*4+y)=M1(:);

end
end



d2logP=zeros(S,16);

P1=1./(B.^2);
P2=1./B;

for x=1:4
for y=1:4

d2logP(:,(x-1)*4+y) = (-1)*( (-1).*P1.*C(:,x).*C(:,y) + P2.*M(:,(x-1)*4+y));

end
end 

d2logP_actual=d2logP(impchoicestate,:);
G=-sum(d2logP_actual,1);

Fisher=reshape(G,4,4);


SE=diag(sqrt(inv(Fisher))/sqrt(nImp));
BETA=[beta_p_sol; beta_x_sol; beta_c_sol; xi_tilde_sol];

[BETA SE]

beta(w,:)=BETA;
standard_errors(w,:)=SE';

end

sum([sign(beta(:,1))  sign(beta(:,1)+1.96*abs(standard_errors(:,1)))  sign(beta(:,1)-1.96*abs(standard_errors(:,1)))   ],2)
sum([sign(beta(:,2))  sign(beta(:,2)+1.96*abs(standard_errors(:,2)))  sign(beta(:,2)-1.96*abs(standard_errors(:,2)))   ],2)
sum([sign(beta(:,3))  sign(beta(:,3)+1.96*abs(standard_errors(:,3)))  sign(beta(:,3)-1.96*abs(standard_errors(:,3)))   ],2)
sum([sign(beta(:,4))  sign(beta(:,4)+1.96*abs(standard_errors(:,4)))  sign(beta(:,4)-1.96*abs(standard_errors(:,4)))   ],2)

cd ../china/estimation/standard_errors
save('standard_errors','standard_errors')


%%%%%%%%%%%%%%%%%%%%%%%%%